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ABSTRACT 


A  fluid-dynamic  computing  method  is  proposed  in  which  the  materials 
are  represented  by  discrete  particles  interacting  with  one  another  by 
means  of  pair  forces.  Details  of  technique,  accuracy,  and  stability  are 
discussed  in  preliminary  form.  Results  are  presented  of  some  simple 
tests  of  the  method,  and  it  is  shown  that  even  though  considerable  de¬ 
velopment  effort  is  yet  required,  the  method  appears  to  have  some  de¬ 
sirable  properties  not  present  in  any  other. 


% 


-3- 


CONTENTS 


Page 


Abstract 

I.  Introduction 

II.  Description  of  the  Method 

A.  Differential  Theory 

B.  Finite  Time  Intervals 

C.  Dissipation 

D.  Neighbors 

E.  Boundary  Conditions 


3 

7 

13 

13 

17 

21 

21 

23 


F.  Fluid-Dynamics  Equivalence,  Choice  of  Force  Functions  24 

III.  PAF  in  One  Space  Dimension  29 

A.  Negative  Internal  Energy  29 

B.  Stability  of  the  PAF  Procedure  30 

C.  Some  Results  of  the  Tests  35 

IV.  PAF  in  Two  Space  Dimensions  40 

V.  Direction  for  Further  Development  43 

Appendix  I  Configurational  Instability  44 

Appendix  II  Angular  Momentum  46 

References  48 


-5- 


I.  INTRODUCTION 


The  perfection  of  high-speed  electronic  computers  has  opened  the 
way  for  solving  a  wide  variety  of  complicated  problems  for  which  no  so¬ 
lutions  previously  were  possible.  Among  these  are  the  problems  of  com¬ 
pressible  fluid  dynamics,  for  which  numerous  techniques  have  been  devel¬ 
oped,  mostly  using  finite-difference  approximations  to  the  appropriate 
partial  differential  equations.  In  many  applications,  particularly  those 
involving  one  space  dimension,  the  results  have  been  very  satisfactory. 
Sufficient  precision  has  been  obtainable  so  that  direct  use  could  be  made 
with  confidence  in  experimental  design  work.  In  situations  depending 
upon  two  space  dimensions,  there  are  likewise  successful  techniques,  but 
at  present  they  all  are  limited  in  applicability.  One  method  may  be 
successful  in  a  situation  in  which  another  would  fail.  In  some  cases, 
problems  have  been  solved  by  using  one  method  for  part  of  the  flow  and 
another  for  the  rest.  This  has  been  quite  successful  for  cases  when  a 
period  of  validity  exists  for  both  methods,  but  has  led  to  frustration 
when  circumstances  would  force  use  of  this  procedure  for  a  problem  in 
which  both  methods  would  give  questionable  results  at  the  time  of  change¬ 
over. 

We  are  here  proposing  a  computing  technique  for  solving  problems  in 
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multidimensional  compressible-fluid  dynamics  that  avoids  some  of  the 
troublesome  features  of  other  existing  methods.  It  is  based  on  a  repre¬ 
sentation  of  the  fluid  by  a  set  of  mass  points  which  are  accelerated  by 
mutual  forces  and  whose  consequent  motions  represent  that  of  a  fluid. 

For  reference,  the  technique  is  referred  to  as  the  "Particle  and  Force" 
(PAF)  method. 

We  illustrate  the  need  for  a  new  computing  method  by  reviewing 
briefly  some  of  the  existing  multidimensional  techniques  which  have  been 
used  at  the  Los  Alamos  Scientific  Laboratory.  This  review  will  also 
serve  as  a  basis  for  demonstrating  some  of  the  properties  of  the  new 
method. 

1 .  Eulerian  finite-difference  methods 

Every  numerical  technique  for  solving  fluid-dynamics  problems  has  as 
its  basis  the  differential  equations  of  motion,  or  equivalently  the  con¬ 
servation  laws  from  which  the  equations  are  derived.  A  principal  differ¬ 
ence  among  techniques  comes  from  the  nature  of  the  coordinate  system  used 
in  the  study.  In  Eulerian  techniques  the  coordinate  system  is  at  rest  re¬ 
lative  to  an  external  observer,  and  fluid  flows  by  the  lines  of  constant 
coordinate.  In  devising  a  finite-difference  approximation  to  the  partial 
differential  equations  of  motion,  one  constructs  a  grid  of  coordinate 
lines  which  divides  space  into  a  set  of  finite  zones  or  cells.  The  fluid 
at  any  given  instant  is  then  described  by  specifying  a  set  of  values  for 
each  cell  of  such  quantities  as  velocity,  temperature,  and  density.  The 
cell-wise  values  are  interpreted  as  some  sort  of  cell-wise  averages. 
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Space  derivatives  in  the  equations  of  motion  are  written  as  the  ratio  of 
intercell  difference  to  cell  size.  In  addition,  the  passage  of  time 
occurs  in  finite  steps  of  interval  &t  each,  so  that  the  time  derivatives 
are  represented  by  the  ratio  of  the  difference  from  end  to  beginning  of 
a  cycle,  to  St.  Thus  the  differential  equations  become  algebraic  equa¬ 
tions  by  which  new  values  for  eacn  cell  (i.e.,  those  at  the  end  of  the 
nth  time  cycle)  can  be  determined  from  the  known  values  at  the  beginning 
of  the  cycle.  Starting  of  a  problem  requires  specification  of  appropriate 
initial  conditions,  and  progress  through  time  requires  boundary  condi¬ 
tions,  translated  into  appropriate  finite-difference  form.  With  proper 
care  in  the  formulation,  and  with  small  enough  space  and  time  intervals, 
the  calculation  may  be  both  stable  and  sufficiently  accurate  ( see  refer¬ 
ence  l).  The  qualification  is  expressed  because  for  some  situations  the 
proper  formulation  is  not  at  present  known,  while  for  others  the  achieve¬ 
ment  of  accuracy  would  be  only  at  the  expense  of  intolerable  amounts  of 
computing  time.  The  most  vexing  difficulties  with  the  Eulerian  methods 
are  of  three  types.  First,  moving  regions  with  structure  which  is  small 
compared  to  the  over-all  dimensions  of  the  system  are  difficult  to  re¬ 
solve  properly.  It  is  wasteful  to  have  fine  zones  throughout  the  entire 
space  through  which  the  small  structure  will  move,  and  it  is  very  diffi¬ 
cult  to  devise  a  good  means  for  creating  and  destroying  zones,  so  that  a 
region  of  fine  resolution  follows  the  structure.  A  second  difficulty 
concerns  the  boundaries  between  materials,  which  must  have  a  special  ex- 

2 

plicit  treatment  or  else  be  subject  to  an  unrealistic  smearing  diffusion.  ’ 
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Third,  there  is  an  inherent  lack  of  Galilean  invariance  and  isotropy. 

One  of  the  overwhelming  advantages  of  the  Eulerian  approach  is  that 
there  is  no  difficulty  in  treating  problems  in  which  there  are  large  dis¬ 
tortions  or  slips  in  the  fluid. 

2.  Lagrangian  finite-difference  methods 

The  approach  is  similar  to  the  Eulerian  described  above,  in  that 
the  space  containing  the  fluid  is  divided  into  finite  zones,  and  the  con¬ 
figuration  advances  through  time  in  finite  steps.  The  fundamental  differ¬ 
ence  is  that  the  Lagrangian  coordinate  system  follows  the  motion  of  the 
fluid. ^  As  a  result,  two  of  the  principal  difficulties  of  the  Eulerian 
methods  are  avoided.  Those  regions  in  which  fine  resolution  is  necessary 
retain  their  fine  zoning  because  the  coordinate  system  moves  right  along. 
Material  interfaces  are  easily  treated  because  coordinate  lines  can  be 
placed  along  them  and  will  forever  follow  them.  The  difficulties  in  the 
Lagrangian  approach  arise  when  the  fluid  develops  internal  slips  or 
large  distortions.  Special  techniques  have  been  developed  to  allow  slip 
lines  to  be  handled;  when  it  is  initially  known,  for  example,  where  the 
slip  will  occur,  then  a  coordinate  line  can  be  placed  there  with  special 
treatment  given.  But  when  the  slip-line  position  is  not  initially  known 
the  matter  is  considerably  complicated.  Distortions  are  even  more  ser¬ 
ious  in  limiting  applicability.  When,  for  example,  initially  rectangu¬ 
lar  cells  have  become  distorted  into  arbitrarily  shaped  quadrilaterals, 
the  initially-appropriate  form  of  the  difference  equations  may  be  quite 
inapplicable,  although,  here  again,  special  techniques  have  been 
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proposed  which  hold  promise  for  solving  the  difficulty  somewhat.  The 
Lagrangian  method  also  has  trouble  with  collapsing  free  surfaces  (as  in 
the  shaped-charge  problem);  some  attempts  to  solve  this  difficulty  are 
also  discussed  in  reference  6. 

3.  Combination  Eulerian- Lagrangian  techniques 

7 

Some  unpublished  work  has  been  done  on  a  technique  in  which  the  co¬ 
ordinate  system  is  Lagrangian  in  one  direction  and  Eulerian  in  another 
( in  two  space  dimensions ) .  Another  alternative  would  be  to  have  a 
Lagrangian  coordinate  system  for  part  of  the  fluid  which  moves  through 
an  Eulerian  system  in  which  the  dynamics  of  the  rest  of  the  fluid  is 
calculated.  Both  approaches  are  nicely  suited  to  some  situations  which 
are  not  well  treated  in  a  homogeneous  coordinate  system  of  either  pure 
type. 

Another  combination  form  which  has  been  extensively  applied  is 

8 

called  the  Particle-in- cell  (PIC)  method.  The  entire  space  occupied  by 
the  fluid  is  covered  by  an  Eulerian  mesh  of  cells,  and  the  fluid  is  re¬ 
presented  by  a  set  of  particles  (essentially  a  Lagrangian  mesh)  which 
move  through  the  cells.  The  finite-difference  equations  are  written  re¬ 
lative  to  the  Eulerian  mesh;  the  particles  are  moved  with  velocities 
determined  from  adjacent  cells.  The  method  has  the  Eulerian  advantage 
of  calculating  distortions  with  ease  and  the  Lagrangian  advantage  of 
giving  information  about  motions  of  fluid  elements,  especially  of  mate¬ 
rial  interfaces  which  are  automatically  handled  well.  PIC-method  disad¬ 
vantages  include  the  Eulerian  difficulties  in  lacking  fine  resolution 
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and  invariance,  together  with  an  additional  disadvantage  of  requiring  the 
storage  of  information  for  both  the  Lagrangian  and  Eulerian  systems. 

4.  Series  methods 

The  dependent  variables  are  expressed  as  infinite  series  in  func¬ 
tions  of  position  with  time-dependent  coefficients.  The  approximation 
comes  in  truncation  of  the  series  so  that  only  a  finite  number  of  coeffi¬ 
cients  remain.  Their  changes  through  time  proceed  in  a  sequence  of  fin¬ 
ite  intervals,  governed  by  equations  obtained  from  substitution  of  the 

series  into  the  appropriate  partial-differential  equations.  Very  little 

9 

work  on  this  method  has  been  done  at  Los  Alamos,  but  Thomas  has  reported 
success  in  application  to  calculations  of  stability  of  plane-parallel  flow. 
A  variety  of  difficulties  concerning  computational  stability,  conserva¬ 
tion,  running  time  on  computer,  and  interpretation  are  forseeable. 

5.  Particle  methods 

The  PAF  method,  which  we  are  proposing  in  this  report,  was  inspired 
by  the  heuristic  studies  of  Pasta  and  Ulam^  on  the  dynamics  of  a  set  of 
particles  with  mutual  forces.  With  forces  dependent  on  separation  only, 
their  calculations  represented  adiabatic  motion  with  no  dissipation. 
Nevertheless,  their  results  on  a  pair  of  Taylor- instability  calculations 
were  quite  encouraging,  and  strongly  suggested  that  a  real  fluid-dynamic 
computing  method  could  be  obtained  from  a  generalization  of  their  ideas. 
Presentation  of  the  resulting  techniques  is  the  main  purpose  of  this 
report . 

11 

Kolsky  recently  has  also  described  a  generalized  particle-like 
method  which  differs  considerably  from  the  one  here  proposed. 
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II.  DESCRIPTION  OF  THE  METHOD 

A.  Differential  Theory 

The  particles  whose  dynamics  we  follow  are  to  represent  a  fluid. 
Insofar  as  we  may  be  aided  in  formulating  the  method  by  referring  to 
classical  particle-dynamic  theory,  we  may  proceed  with  that  guidance.  We 
know,  however,  that  at  some  point  a  divergence  will  be  necessary  so  that 
the  dissipative  effects  in  a  real  fluid  can  be  represented.  Our  particles 
are  not  molecules  whose  internal  energy  is  carried  by  velocity  fluctua¬ 
tions;  indeed,  we  expect  that  the  velocity  of  a  particle  is  to  represent 
the  mean  velocity  of  the  finite  mass  of  fluid  it  represents.  The  macro¬ 
scopic  kinetic  energy  of  the  fluid  is  exactly  the  kinetic  energy  of  all 
the  particles.  The  internal  energy  is  represented  by  an  additional  var¬ 
iable.  If  this  latter  is  expressed  as  a  function  of  particle  positions 
only,  then  only  adiabatic  motion  can  be  represented.  Compression  and 
subsequent  expansion  can  return  the  set  of  particles  to  exactly  their 
initial  configuration  with  no  dissipation.  Thus  a  special  prescription 
is  needed  to  describe  variations  of  particle  internal  energy. 

We  start  by  considering  the  dynamics  of  a  set  of  particles  de¬ 
scribed  by  the  following  nomenclature. 
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i,j  h  indices  describing  the  particle  number 

m.  s  mass  of  particle 
J 

"r.  a  space  coordinate  of  particle  //j 
J 

u.  =  velocity  of  particle  #o 
J 


F.  .  =  force  exerted  by  particle  i  onto  particle  j 
i  J 


— >  — >  — > 

r .  .  =  r .  -  r 
i<3  J  1 


r.  .  = 
U 


r.  . 
U 


=  ?.  ,/r.  .  (a  unit  vector  pointing  from  particle  #i 
ij  iO  i0  particle 


M  =  m  .u .  =  momentum  of  particle 
0  0  J 

K  5-b  u  .u.  s  kinetic  energy  of  particle 
0  2  j  j  j 


J 


j 


=  internal  energy  of  particle 


#0 


Additional  nomenclature  will  be  introduced  as  required  by  developments 
of  the  theory. 

We  commence  by  assuming  that  the  particles  are  governed  by  the 
equations  of  motion 


i 


0) 


(2) 


The  summation  over  i ,  modified  by  the  presence  of  does  not  include 
the  term  i  =  j,  and  is  further  restricted  to  include  only  certain  neigh¬ 
bors  of  j  as  discussed  further  ahead.  Summation  without  *  includes  all 
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particles  in  the  system 


Next,  we  assume  that  the  force  function  can  be  divided  into  two 


parts 


s. .f.  .  + 
U  10 


10 


(3) 


where 


f .  . 
10 


Ji 


The  first  term  in  the  force  is  associated  in  form  with  the  equation 
of  state  of  the  fluid;  the  second  term  is  introduced  to  achieve  dissipa¬ 
tion  in  the  same  manner  as  the  "artificial  viscosity"  of  von  Neumann  and 
1  P 

Richtmyer,  or  for  the  purpose  of  including  real  viscous  effects. 

The  correspondence  with  fluid  mechanics  comes  through  an  examination 
of  the  conservation  laws  in  forms  appropriate  to  the  nature  of  the  con¬ 
tinuum  to  be  represented. 

1.  Conservation  of  mass.  This  is  automatic.  Each  particle  has 

constant  mass,  m.,  so  that  the  total  does  not  change  with  time.  Like- 
<3 

wise  the  change  of  mass  in  any  fixed  volume  exactly  equals  the  amount 
flowing  over  the  bounding  surface. 

2.  Conservation  of  momentum.  To  satisfy  this  requirement,  the  re¬ 
striction  is  the  same  as  in  classical  particle  dynamics,  namely 

— >  — > 

F. .  =  -  F  .  Proof  of  this  is  demonstrated  as  follows.  Consider  the  mo- 
ij  ji 

mentum  change  rate  of  a  particular  subset  of  all  the  particles 


j( subset)  j( subset)  i 
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We  break  the  sum  over  i  into  two  parts  and  write 


i-  r  l  fu-  I  I  n, 

i(inside)  j( subset)  i(outside)  j( subset) 

where  "inside"  and  "outside"  refer  to  inclusion  or  exclusion  from  the 
subset.  In  the  first  double  sum,  each  pair  of  particles  enters  twice, 
so  that  the  total  contribution  to  the  sum  from  a  particular  pair  is 
if  +  F*  .  Since  there  must  be  no  contribution  to  momentum  change  from 

ij  ji 

particles  within  the  subset,  the  sum  of  the  two  terms  must  vanish.  The 
second  double  sum  does  not  thereby  vanish,  since  each  pair  enters  only 
once.  Thus  with  F4  .  =  -  F_.  ^ ,  the  momentum  change  of  any  subset  of  par- 
tides  arises  only  through  external  forces,  as  required. 

The  restriction  also  means  that 

f. .  e  f .. 
ij  Ji 

_  h 
gij  -  '  g 

3.  Conservation  of  energy.  Here  we  must  make  a  break  from  the 
usual  procedure  in  particle  mechanics.  In  devising  a  reasonable  approach 
we  will  thereby  be  able  to  establish  some  of  the  crucial  parts  of  the 
technique.  The  basis  for  the  energy  discussion  is  that  the  rate  of 
change  of  energy  of  a  particle  should  be  given  by  the  rate  that  the 
other  particles  do  work  on  it.  This  work  rate  is  in  turn  given  by  the 
product  of  force  by  velocity.  (To  be  properly  symmetric,  the  velocity 
through  which  the  work  flux  is  carried  from  one  particle  to  another 
must  be  the  mean  value  of  the  two  velocities.)  Thus  we  write 
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(5) 


_d_ 

dt 


* 

(K,i +  V  ■  I 

i 


F.  .* 
ij 


It  follows  that  the  total  energy  of  an  isolated  system  is  conserved 

* 


dE 

dt 


_d_ 

dt 


I(Kj +  V 


HI 


— >  — >  — ) 

F.  .*(u.  +  u.)  =  0 
ij  i  J 


where  equality  to  zero  follows  from  separate  vanishing  of  the  sum  of  con¬ 
tributions  from  each  pair.  Likewise  the  energy  of  any  subset  of  particles 
changes  only  through  work  done  on  them  by  external  particles. 

Now,  we  already  know  that 


dX, 

dt 


__a  =  -  * 


du. 


m.u.«  ,, 

J  J  dt 


-tf.-V?.. 

j  L>  ij 


Combination  of  this  with  (5)  can  be  arranged  to  give 


dJ. 

dt 


1  y  f  .(u . 

2  L  ij  V  i 


V 


(6) 


or 


dJ. 

_ i 

dt 


'  -  5  Z  *  5  X  V(?i  ‘ 


(7) 


B.  Finite  Time  Intervals 

In  practice,  the  numerical  computations  must  proceed  through  a 
sequence  of  finite  time  advancements,  whose  steps  are  of  duration  5t. 
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This  is  accomplished  through  a  replacement  of  (l),  (2),  and  (6)  by 


->  n+1  -» n  * 

u .  -  u .  V1 


m 5sL— 
3 

6t 

-  -  L  *±3 

i 

n+1 
r . 

J_ 

n 

-  r . 

3.  .  . 

n+1 

—  n 

St 

—  U-  . 

3 

n+1 

u  . 

_J _ 

-  jn 

,  A  „ 

.1  V*?« 

5t 


>  -  <  u .  > 
3 


where 


.  — i  /— »n  -+  n+1  \ 
<  U  >  =  5  (u  +u  ) 


(8) 

(9) 


(10) 


This  shows  how  the  variables  for  the  beginning  of  cycle  #n+1  are  ob¬ 
tained  algebraicly  from  those  at  the  beginning  of  cycle  #n.  The  choice 
of  time-centering  of  the  equations  is  justified  as  follows: 

Eq.  (8)  —  At  the  beginning  of  the  calculations  for  the  advance¬ 
ment  through  a  cycle,  the  only  information  available  for  the  force  cal¬ 
culation  is  that  which  pertains  to  the  beginning  of  cycle  #n.  The  force 
is  thus  labeled  with  index  n. 

Eq.  (9)  —  After  calculation  of  the  new  velocity  by  (8),  there  is 
some  arbitrariness  as  to  what  velocity  should  be  used  to  determine  the 
new  coordinate.  To  show  that  the  newly  calculated  velocity  is  preferred 
over  that  which  the  particle  had  at  the  beginning  of  the  cycle,  we  appeal 
to  a  simple  example  of  stability  properties.  Suppose  that,  in  one  dimen¬ 
sion,  a  particle  is  subjected  to  an  external  potential  which  is  a  function 
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of  the  particle  position  only.  Then 


dx 

dt 


are  the  appropriate  equations  for  small  oscillations  about  a  potential 
minimum  centered  at  x  =  0.  The  solution  can  be  written 
itD  t 

u  =  u  e  ° 
o 

t 

x  =  x  e  ° 
o 


In  finite -difference  form,  analogous  to  (8)  and  (9),  we  write 


n+1  n  2  „  n 

u  -  u  =  -  CD  St  x 
o 

n+1  n  n+e 
x  -  x  =  u  5t 


where  if  e  -  1  we  have  the  proposed  procedure,  and  if  e  =  0  we  have  the 
alternative  whereby  the  particle  is  moved  with  its  beginning-of-cycle 
velocity.  We  try  the  solution 

n  icjon&t 

u  =  u  e 
o 

xn  =  x  eicon&t 
o 

and  find  the  condition  for  solution  to  be 
(r  -  1  )2  =  -  re(cuo6t)2 
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where  r  =  ei03&t.  Instability  of  the  difference  approximation  is  indi¬ 
cated  by  the  presence  in  co  of  a  negative  imaginary  part,  or  equivalently 
by  a  magnitude  of  r  which  exceeds  unity.  For  e  =  0,  r  =  1  ±  ic00&t  and 
|r |  =  +Cn25tH  ,  indicating  unconditional  instability.  For  e  =  1, 

however,  r  -  (<e^)2  */  («»)4  7  “  “o6t  »  2’  ““  " 

is  real  and  one  of  the  solutions  has  magnitude  exceeding  unity.  If 
ffi()&t  <  2,  then  |r|  ■  1,  indicating  the  achievement  of  stability  for  suf¬ 
ficiently  small  values  of  6t.  (The  time  interval  per  cycle  must  be  less 

than  n”1  times  the  period  of  oscillation.) 

Eg.  (10)  -  The  right  side  contains  the  average  of  the  old  and  new 

velocities,  which  combination  is  introduced  to  assure  rigorous  energy 
conservation  in  the  time-difference  form  of  the  equations.  (Mass  and 
momentum  have  likewise  been  conserved;  proof  of  this  is  the  same  as  for 
the  differential  equations.)  To  demonstrate  energy  conservation  we  start 

from  the  identity 


2 


fer)2  -  (*" 


_»  /_» n+1 

=  <  U.  >  (  u.  -  U. 

J  \  0  j. 


Thus,  from  (8),  the  change  in  kinetic  energy  of  a  particle 


* 

V  -» 


Kj+1  -  K/  =  &t  L  <  ^ 


Combination  of  this  with  (10)  gives  for  the  change  in  total  particle 


energy,  E  , 
J 


-,n+1 


-EJn'5l?ir(<“i>  +  <?J>) 


(ID 
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Since  F.  .  =  -  F.  this  result  shows  that  the  energy  transferred  from 
ij  ij 

particle  # i  to  particle  j  is  equal  in  magnitude  but  opposite  in  sign 
from  that  transferred  from  to  #i,  thus  proving  the  contention  of 
conservation. 


C.  Dissipation 

There  is  at  least  one  other  form  alternative  to  (10)  which  could 
be  considered  for  the  internal  energy  calculation.  The  total  energy 
difference  could  be  calculated  in  a  form  analogous  to  ( 11 ) ,  but  with  any 
time  centering  of  the  right  side.  The  result  would  still  be  conserva¬ 
tive  as  long  as  the  proper  reciprocal  symmetries  were  preserved.  From 
the  new  total  energy  of  the  cell,  the  new  kinetic  energy  could  then  be 
subtracted  giving  the  new  internal  energy. 

The  reason  for  the  specific  choice  of  the  form  (10),  however,  fol¬ 
lows  from  the  requirement  of  monotonic  dissipation.  With  (lO),  proper 
choice  of  the  forces  can  always  result  (at  least  to  lowest  order  in 
5t )  in  increasing  entropy,  while  in  most  of  the  alternative  forms  there 
can  be  circumstances  leading  to  significant  decreases  in  entropy. 

D.  Neighbors 

Two  kinds  of  particle-wise  sums  have  been  introduced  so  far.  The 

unmodified  summation  is  over  all  particles  of  the  system,  while  a  summa- 

*  -> 

tion  modified  by  superscript  *  means  a  neighbor  sum.  Thus  L  F_  de¬ 
notes  a  sum  over  the  neighbors,  i,  of  particle  Such  a  modified  sum¬ 

mation  is  required  in  order  to  achieve  a  proper  fluid-dynamic  represen¬ 
tation,  in  which  elements  of  fluid  do,  indeed,  interact  only  with 
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those  adjacent.  There  are  several  methods  by  which  the  neighbors  of  a 
particle  could  be  chosen,  so  that  some  criteria  must  be  evolved  for  de¬ 
ciding  among  them.  Important  among  these  criteria  is  the  requirement 
of  reciprocity:  for  rigorous  conservation  of  momentum  and  energy  it  is 
necessary  that  particle  |i  be  a  neighbor  of  particle  $j  if  and  only  if 
#j  is  a  neighbor  of  #i.  This  is  satisfied  by  the  neighbor-choosing 
method  in  which  all  those  particles  lying  closer  than  some  given  dis¬ 
tance  to  #j  are  considered  to  be  neighbors  of  #j.  The  distinguishing 
distance  must  be  the  same  for  all  particles,  with  the  result  that  in  re¬ 
gions  of  high  compression  each  particle  will  have  many  neighbors,  in 
violation  of  the  requirement  of  interaction  only  with  adjacent  elements. 
We  have  therefore  used  a  somewhat  different  means  for  finding  neighbors. 

The  neighbor-choosing  procedure  we  have  used  involves  two  steps. 
First,  for  each  particle  we  find  the  N  other  particles  lying  closest  to 
it.  (in  Cartesian  coordinates,  N  is  twice  the  number  of  space  dimen¬ 
sions.)  Since  the  reciprocity  condition  is  not  automatically  satisfied 
by  this,  we  proceed  in  the  second  step  to  add  to  the  neighbors  of  each 
particle  just  those  necessary  to  give  complete  reciprocity.  As  a  re¬ 
sult,  each  particle  can  have  more  neighbors  than  the  ideal  number,  N, 
but  we  do  not  therefore  modify  the  nature  of  the  interparticle  forces 
or  any  of  the  procedures  so  far  described.  An  argument  justifying  this, 
and  a  demonstration  of  its  validity,  are  given  in  Part  III. 

The  following  table  shows  how  many  neighbors  a  particle  would 
have  under  various  circumstances,  in  Cartesian  coordinates. 
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Number 

of 

Dimensions 

72 

s  Number  of  Neighbors 

N 

Ordinary 

Total 

Maximum  Possible 
Total 

Usual  Extreme 
Total 

1 

2 

~  2 

4 

3 

2 

k 

~  5 

9(?) 

8 

3 

6 

~  8 

15(?) 

12 

The  numbers  are  not  meant  to  have  any  absolute  significance,  only  to  indi¬ 
cate  approximately  the  magnitudes  to  be  expected.  The  numbers  for  three 
dimensions  were  guessed  as  extrapolations  from  one-  and  two-dimensional 
observations. 

A  combination  of  the  above  two  procedures  also  might  be  appropriate. 
The  search  for  N  nearest  neighbors  could  be  conducted  only  among  those 
which  lie  closer  than  a  prescribed  distance.  This  would  mean  that  a  small 
number  of  particles  far  from  the  rest  could  actually  become  completely  de¬ 
tached,  and  no  longer  influence  the  rest. 

E.  Boundary  Conditions 

The  boundary  adjacent  to  a  vacuum  is  simply  a  point,  line,  or  sur¬ 
face  separating  a  region  with  particles  from  one  without  any.  One  might 
speculate  that  no  special  consideration  would  be  necessary  for  obtaining 
properly  the  surface  motion,  and  in  Part  III  it  is  shown  that  such  is 
indeed  the  case. 

A  rigid  wall  can  be  produced  by  an  image  method  in  which  the  neigh¬ 
bors  of  a  particle  near  the  wall  include  the  images  of  itself  and  of  its 
neighbors.  Each  image  particle  has  the  reflected  properties  of  the  cor¬ 
responding  one  in  the  system.  This  is  easily  accomplished  if  the  wall 
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is  straight  (or  flat),  but  becomes  more  complicated  if  the  wall  is  curved. 
In  this  latter  case  one  might,  for  simplicity,  be  content  to  allow  a 
given  particle  to  be  influenced  by  its  own  image  only  (or  images,  if  the 
wall  forms  a  small  pocket). 

So  far  our  studies  have  been  limited  to  examples  with  only  the  sim¬ 
plest  boundary  conditions;  much  thought  and  experimentation  will  be  re¬ 
quired  to  perfect  general  techniques. 

F.  Fluid-Dynamics  Equivalence.  Choice  of  Force  Functions 

The  most  difficult,  but  most  important,  question  to  be  answered  is: 
How  well  will  results  from  the  proposed  computing  method  represent  the 
solution  of  a  fluid-dynamics  problem?  At  this  time  only  a  partial  answer 
can  be  given,  and  most  of  the  rest  of  this  report  will  show  that  which  so 
far  has  been  learned. 

Up  to  this  point,  no  specification  has  been  made  of  the  force  func¬ 
tion,  except  as  implied  by  the  dependence  of  the  s*  .f .  .  part  upon  the 
mutual  distance,  the  internal  energies,  and  an  ordering  variable  (such 
that  the  force  vanishes  between  two  particles  if  they  are  not  neighbors), 
together  with  the  statement  that  the  g^  part  is  included  only  to  produce 
dissipation.  Already  it  follows  that  goodness  of  the  method  will  depend 
upon  some  statistical  properties  of  the  particle  behavior.  An  ideal 
particle-and-force  treatment  might  base  the  force  between  two  particles 
upon  particle-wise  volumes  and  interparticle  area.  (Both  of  these  depend 
in  complicated  fashion  upon  the  positions  of  all  the  other  surrounding 
particles,  and  would  therefore  require  very  large  amounts  of  computing 
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time. ^  In  the  PAF  technique,  however,  the  force  between  two  particles  is 
completely  independent  of  the  other  surrounding  particles,  with  the  im¬ 
portant  exception  that  the  two  particles  must  be  neighbors. 

The  s*.  .f  .  part  of  the  force  function  has  been  chosen  in  such  a  way 

U  id 

as  to  represent  the  equation  of  state  in  a  rectangular  array  of  particles. 
Consider  the  case  of  two  dimensions,  with  the  particles  all  having  the 
same  mass,  m.  For  a  polytropic  gas,  for  which  the  equation  of  state  is 
p  =  (y  -  l)pl  (in  which  p,  p  and  I  are  respectively  pressure,  density,  and 
specific  internal  energy),  the  force  between  two  particles 
in  a  horizontal  direction  is  thus 


The  same  form  results  for  the  vertical  direction,  and  (12)  is  thus  the 
force  function  (actually  appropriate  to  use  in  one,  two,  or  three  dimen¬ 
sions)  . 

More  difficult  to  derive  is  the  form  which  should  statistically  re¬ 
present  a  more  general  equation  of  state.  For  p  =  p(p,l),  the  recipe  cor¬ 
responding  to  the  polytropic  gas  treatment  would  lead  to 
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horizontal: 

vertical: 


The  result  cannot,  therefore,  be  directly  related  to  the  requirement 
which  states,  in  effect,  that  the  horizontal  force  should  not  depend 
upon  b,  while  the  vertical  force  should  not  depend  upon  a.  (if,  however, 
the  motion  is  purely  horizontal,  then  b  is  constant  and  the  requirement 
is  satisfied. ) 

Actually,  a  rectangular  array  is  configurationally  unstable,  as 

14 

has  been  pointed  out  by  Birkhoff  and  Lynch,  and  the  particles  will 
most  usually  be  distributed  "randomly."  (A  discussion  of  configuration 
instability  is  given  in  Appendix  I.)  In  addition,  the  randomness  of 
spacing  will  probably  be  maintained  at  least  at  a  rate  comparable  to  the 
rate  of  the  over-all  dynamics.  Thus  the  area  between  particles  should, 
statistically,  be  proportional  to  the  mean  interparticle  distance,  while 
the  volume  used  to  calculate  the  density  should  be  proportional  to  the 
square  thereof.  This  leads,  therefore,  to  the  proposal  for  general  equa¬ 
tion  of  state 


_  a-i 

f .  .  =  r. .  p 
ij  U 


m 


J.  +  or. 

i _ o 


L(r.  .) 
10 


,cr 


2m 


(13) 


to  be  appropriate  for  ol  space  dimensions.  (The  constants  of  proportion¬ 
ality  are  chosen  to  give  agreement  with  the  polytropic  equation  of  state, 
and  the  density  formula  for  a  square  array.) 


1 5 

Incidentally,  Beyer  has  pointed  out  that,  at  least  with  a 
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rectangular  array  of  particles,  an  infinite  line  at  any  angle  through 
them  will  have  the  proper  force  acting  across  it  when  the  fluid  is  a 
polytropic  gas.  It  would  be  useful  to  have  a  generalization  of  his  ob¬ 
servation  to  the  determination  of  the  probability  of  finding  a  specified 
force  across  a  finite  line  segment  through  a  random  distribution  of  par¬ 
ticles  with  generalized  force  function. 

The  g  .  forces  are  much  more  easily  specified  in  general.  Indeed, 

1 J 

we  have  chosen 


— >  f  — >  —4  \ 

Hi  “  (ui  -  V 


(14) 


for  a  region  of  compression  and  g. .  =  0  for  a  region  of  expansion;  o>  is 

a  constant  with  dimensions  of  reciprocal  time.  (An  alternative  form  is 

discussed  in  Appendix  II.)  The  analogy  to  viscosity  is  easily  demon- 

»■) 

strated  and  the  contribution  that  the  g^j  force  makes  to  stability  is 
discussed  in  Part  III.  The  dynamic  effect  can  be  seen  as  follows.  Con¬ 
sider  a  group  of  particles  for  which  the  f. .  forces  vanish.  Then  (l) 
becomes 

du.  * 

m  =  mcu  Z  (tu  -  u*. ) 
dt  i  i  j 


where  xx. 1  =  —  Y,  u^  is  the  average  velocity  of  the/Z  neighbors.  If  this 


J  'Til 
cor 

will  change  velocity  according  to 


average  is  constant  in  time  then  particle  jfj ,  starting  from  u.  at  t  =  0, 

JO 


u .  =  u . 1  +  (u  .  -  u .  * )  e 

J  J  jo  J 
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The  velocity  will  decay  to  the  mean  velocity  of  its  neighbors,  with 
a  decay  time  of  (?2a>)  1 . 

Thus,  the  crudeness  of  our  proposed  method  is  now  laid  bare.  If  it 
can  be  made  to  work  as  proposed,  or  with  relatively  minor  variations,  then 
it  will  have  tremendous  advantages  of  generality  without  requiring  the 
great  computing  time  of  most  proposals  of  this  generality.  The  proof  of 
applicability  must  come  through  considerable  analytical  study  of  the  sta¬ 
tistics  of  such  a  system,  together  with  an  analysis  of  numerous  computing 
experiments.  These  projects  are  in  progress  in  Group  T-3  of  the  Los 
Alamos  Scientific  Laboratory,  and  contributions  or  suggestions  will  be  en¬ 
thusiastically  welcomed. 
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III.  PAF  IN  ONE  SPACE  DIMENSION 

The  one-dimensional  form  of  PAF  differs  little  from  the  usual  tech¬ 
niques  of  one -dimensional  Lagrangian  finite  difference  methods.  The  prin¬ 
cipal  difference  for  which  concern  was  felt  is  the  fact  that  a  PAF  par¬ 
ticle  can  have  either  one  or  two  neighbors  on  each  side,  as  a  result  of 
the  reciprocity  requirement.  Thus,  for  example,  as  a  compression  wave 
moving  to  the  right  approaches  a  particle,  there  will  come  a  time  when 
that  particle  will  have  two  neighbors  to  the  left  and  one  to  the  right. 
Briefly,  its  rightward  acceleration  will  be  too  great.  Soon,  however, 
the  situation  is  reversed  —  after  the  front  of  the  compression  wave  has 
passed,  the  particle  will  have  for  another  brief  time  two  neighbors  to 
the  right  and  only  one  to  the  left.  As  a  result,  the  effects  of  previous 
over-acceleration  can  be  expected  to  be  removed;  the  computations  dis¬ 
cussed  below  show  that  the  expectation  is  indeed  realized. 

The  one- dimensional  calculations  were  performed  for  the  additional 
purposes  of  studying  the  behavior  of  internal  energy  and  certain  matters 
concerning  stability. 

A.  Negative  Internal  Energy 

With  the  neighbor-search  procedure  which  gives  each  particle  at 
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least  N  neighbors,  even  when  it  has  flown  off  from  the  rest  of  the  system, 
it  is  possible  to  obtain  negative  values  for  the  specific  internal  energy. 
This  can  be  seen  by  considering  the  motion  of  a  particle  which  is  at  a 
distance  x  from  its  two  neighbors  (both  of  which  are  on  the  same  side  of 
the  particle).  If  the  particle  has  internal  energy  J,  and  its  neighbors 
have  internal  energy  JQ  (assumed  constant),  then,  for  J  small, 

dJ  (7-0Jou  (7-OJo1dx 

dt  2x  "  2  x  dt 

(y  -  DJ 

J  ~  const - g -  In  x 

and  since  x  continues  to  increase  (since  the  force  does  not  become  attrac¬ 
tive  until  somewhat  after  J  =  0),  it  is  seen  that  J  could  become  negative. 
This  situation  could  probably  be  remedied  by  the  combination  neighbor¬ 
searching  procedure  described  at  the  end  of  Part  II-D. 

B.  Stability  of  the  PAF  Procedure 

Consider  a  one -dimensional  region  which  is  only  slightly  perturbed 
from  a  constant  state.  For  simplicity,  we  use  a  polytropic-gas  equation 
of  state,  but  the  results  are  more  generally  applicable.  The  equations 
of  motion  are 
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n+1  n  n+1  ~ , 

xj  -  xj  =  u3  &t 

The  values  of  £+  are  zero  and  one,  for  expansion  and  compression  respec¬ 
tively.  Into  these  equations  we  introduce  the  linearizations 

x.  -  x.  s  6x(l  +  e  i) 

0  0-1  0-2 

j  a  J(1  +  2 a  ) 

J  J 

and  drop  higher  than  first  order  terms  in  u.,  e  iand  a  ;  J  and  6x  are 

J  J~z[  J 

constants.  The  resulting  linearized  equations  of  motion,  plus  energy 
equation,  are 


m  f  n+1  n^  (y  -  1  )J  ( _n  „n  ,  jn  en  . 

rr  u.  -  U .  J  =  r - —  l  0  .  -  d.  ,  t  +  e  .i  ~  €  i  JL 

5t  \  o  0/  \  0-1  j+1  0+2  0 " 2 


.  (  n  n 

'iMVi  -  ui 


u  (  n  n 

j+i  j. 


6t  /  n+1  n+1  A 

Vh  •  Viy 


-  1 )  &t  /  n 
8  5x  VUj-1 


n  n+1  ' 

•  uj+i  •  uo+iy 


We  assume  that  the  time  average  of  £+  is  \  and  put  this  value  in.  To 
analyze  the  stability  of  these  equations  we  substitute  the  trial  func¬ 


tions 


n  ikj  n 

u .  s  u  e  r 
J  o 


n  ikj  n 

e  .  ==  e  e  r 
J  o 
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n  ikj  n 

a  =  o  e  r 

J  o 


and  find  them  to  he  a  solution  provided  that  r  satisfies  the  equation 

2 

(r-l)[r-1+|(l  -  cos  k)]  +  y-  [4r  sin2  ^  -  (l  +  r)  sin2k]  =  0 


in  which 


|  =  co&t 
cSt 


P  = 


6x 


and  the  sound  speed  is  c  s  \fy(y  -  1 ) J/m  .  The  equations  are  not  stable 
for  those  circumstances  in  which  the  magnitude  of  r  exceeds  unity.  Con¬ 
sider  first  the  case  in  which  the  roots  of  r  are  real.  At  r  =  1  we  find 
that  except  for  certain  peculiar  values  of  k,  it  is  necessary  that  p  =  0, 
while  for  p  >  0,  r  <  1,  and  no  restriction  therefore  results.  At  r  =  -1, 


\^_ _ 2 

^  +  7  1  -  cos  k 


If  |  +  h2/7  exceeds  the  right  side,  then  r  <  -  1 .  The  most  restrictive 
case  comes  from  the  smallest  right  side,  and  can  be  expressed  in  the  sta¬ 


bility  requirement 
2 

i  +  —  <  l 
7 

Finally,  we  must  examine  the  complex  roots  of  r,  which  can  be 
shown  to  have  magnitude 


(15) 


|r  |  =  y  1  -  6(1  -  cos  k)  +  y-2-~-  —  V2  sin2 
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so  that,  at  |r|  =  1, 


£  =  ~2~~  ^0  +  cos  k) 

In  this  case,  if  £  is  less  than  the  right  side,  then  ]r|  >  1.  As  a  func¬ 
tion  of  k,  the  right  side  has  maximum  value  at  cos  k  =  1 .  We  thus  obtain 
the  additional  restriction  for  stability 

(16) 

7 


In  more  familiar  nomenclature,  the  two  stability  conditions  can  be  sum¬ 
marized 


CD  5t 


6t  < 


(17) 


The  result  may  be  compared  to  that  for  Example  3,  p.  18  of  Ref.  1,  in 
which  a  similar  stability  analysis  is  made.  Finally,  with 


11 


_  cSt 
~  5x 

cuSx 

c 


being  dimensionless  measures  of  5t  and  cd,  we  may  write  the  stability 
conditions 


1  2 

W  +  y  ^  <  1 

M  r]<  0 


08) 
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As  an  example.  Fig.  1  shows  a  plot  of  the  stable  region  for  7  =  5/3- 


FIG  .  1 

The  points  were  tested  by  machine  calculations;  those  marked  •  were 
found  to  be  stable  while  those  marked  x  were  found  to  be  unstable. 

The  results  of  this  stability  analysis  are  similar  to  those  which 
would  result  for  a  much  more  general  class  of  equations  with  both  wave- 
like  and  diffusion-like  properties.  In  general,  the  conditions  for  sta¬ 
bility  can  be  violated  in  either  of  two  ways  (corresponding  to  leaving 
the  stability  region  either  to  the  right  or  to  the  left,  in  Fig.  1).  In 
one  case,  there  is  too  much  diffusion  for  the  given  sizes  of  5t  and  &x. 

A  perturbation  to  an  otherwise  smooth  profile  is  over-corrected,  and  the 
result  is  oscillation  with  increasing  amplitude.  Such  instability  is 
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easily  observed,  and  is  the  nature  of  all  the  unstable  points  marked  in 
Fig.  1 .  The  other  type  of  instability  corresponds  to  too  little  diffu¬ 
sion.  Since  the  equations  without  the  diffusion  properties,  and  with 
6x  -»  0,  8t  -> 0,  have  solutions  with  constant  amplitude,  the  truncation 
resulting  from  finite  differences  can  easily  be  expected  to  add  slowly- 
growing  influences.  Such  is  the  case  with  the  unstable  region  to  the 
left  in  Fig.  1;  and  while  two  attempts  were  made  to  observe  the  insta¬ 
bility  in  that  region,  the  rate  of  growth  was  so  small  as  to  leave  doubt 
concerning  confirmation  of  instability.  Actually  these  results  are  not 
unreasonable  since  the  equation  for  r,  giving  the  rate  of  growth,  shows 
that  the  examples  tried  (both  chosen  for  large  growth  rate)  correspond 
to  two  too  small  rates  to  be  observed  with  the  computing  code  used. 

C.  Some  Results  of  the  Tests 

1 .  The  rarefaction  wave.  Gas  adjacent  to  a  wall  at  x  =  0  is 
initially  at  rest.  Beyond  the  gas  to  the  right  is  vacuum.  Initial  data 


for  the  PAF  calculation  are 

Number  of  particles  =  30 

Internal  energy  per  particle  =  1 .0 

Specific  heat  ratio,  y  =2.0 

Mass  per  particle  =  1 .0 

6x  =1.0 

6t  =  0.1 

oo  =1.0 
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The  initial  sound  speed  is  thus  1.4l4  and  the  escape  speed  of  the 
free  surface  is  2.8l8.  The  configuration  of  particles  together  with  their 
internal  energy,  velocity,  and  specific  volume  (i.e.,  8x)  are  shown  as 
functions  of  position  for  time  t  =  20  in  Fig.  2.  The  datum  points  are 
from  the  machine  calculation,  while  the  solid  lines  show  the  true  solution. 


FIG.  2 

2.  The  shock  wave.  Gas  adjacent  to  a  wall  at  x  =  0  is  initially 
cold  and  moving  towards  the  wall.  Initial  data  for  the  PAF  calculation 


Number  of  particles 

=  51 

&x 

=  i  .0 

Velocity  of  particles 

=  -  2.0 

6t 

=  0.1 

Specific  heat  ratio,  7 

«  2.0 

(JD 

=  1 .0 

Mass  per  particle 

=  1 .0 
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are 


The  configuration  of  particles,  together  with  their  internal  energy, 
are  shown  for  time  t  =  1 5  in  Fig.  Datum  points  are  from  the  machine 
calculation,  while  the  solid  line  shows  the  true  solution.  Note  the  usual 
Lagrangian  difficulty  of  overproduction  of  internal  energy  at  the  wall. 


The  results  from  the  first  shock  calculation  are  also  shown  as  a 


function  of  time  in  Fig.  5.  The  graph  shows  the  total  energy  of  the  sys¬ 
tem  together  with  the  negative  of  the  total  momentum.  Datum  points  are 
from  the  calculation  while  solid  lines  show  the  true  solution. 


Finally,  in  Fig.  6  is  shown  as  a  function  of  time  the  difference 
between  calculated  and  true  histories  of  internal  energy  for  various  values 
of  co. 
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TIME 


FIG.  6 

It  is  apparent  that  the  overproduction  arises  at  early  times  when 
the  shape  of  the  smeared  shock  is  adjusting  itself.  Once  this  is  com¬ 
pleted,  subsequent  internal  energy  production  proceeds  at  very  nearly  the 
proper  rate. 
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IV.  PAF  IN  TWO  SPACE  DIMENSIONS 


The  only  two-dimensional  tests  so  far  performed  have  been  of  a  re¬ 
latively  simple  nature,  for  the  purpose  of  revealing  whether  or  not  some 
basic  difficuty  exists.  The  test  results  have  been  encouraging,  however, 
and  we  are  therefore  developing  a  much  more  elaborate  computing  code 
which  will  allow  the  interaction  of  several  materials  with  several  boun¬ 
daries. 

In  the  simple  tests,  a  region  of  gas  was  allowed  to  interact  with 
one  rigid  wall.  In  the  first  calculation,  the  gas  region  was  rectangular, 
with  one  side  initially  along  the  wall,  and  the  gas  was  cold  but  possessed 
a  velocity  towards  the  wall.  Initial  conditions  were  much  like  those  in 
the  one-dimensional  shock  problem. 


Number  of  particle  rows 

(parallel  to  wall)  =  8 

Number  of  particle  columns 

(normal  to  wall)  =  15 

Internal  energy  per 

particle  =  0 

Velocity  of  particles 

(parallel  to  wall)  =  0 


Velocity  of  particles 

(normal  to  wall)  =  -2.0 

Mass  per  particle  =  1 .0 

Specific  heat  ratio, 7  =  2.0 

&x  =  5y  =1.0 

5t  =0.1 

co  =1.0 
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The  problem  is  one  dimensional  along  the  center  line  until  signals  from 
the  side  arrive.  The  two-dimensional  effects  come  from  lateral  splash¬ 
ing  of  the  rectangle  of  gas.  While  exact  comparison  solutions  for  the 

two-dimensional  effects  were  not  obtained,  it  was  possible  to  approximate 

16 

the  effects  using  simple  elements  of  shock  and  rarefaction  theory. 

Thus,  Fig.  7  shows  the  time  history  of  the  normal  momentum  of  the  gas. 

The  solid  line  is  the  theoretical  one-dimensional  history  (that  is,  the 
history  without  splash  effects),  while  the  dashed  line  shows  the  lowest 
order  solution  with  splash  effects  included.  Datum  points  are  from  the 
calculation. 


FIG.  7 
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Considering  the  crudeness  of  the  computational  resolution,  the  agreement 
appears  quite  good.  Figure  8  shows  the  tangential  momentum  of  one  side 
of  the  gas. 


FIG.  8 


Here  the  effect  is  entirely  two-dimensional,  associated  with  the  splash. 
The  theoretical  curve  (solid  line)  is  based  upon  the  lowest  order  effects. 
Later  times  were  not  shown  because  the  configurational  instability  ( see 
Appendix  I)  mentioned  earlier  made  determination  of  the  computed  momen¬ 
tum  difficult.  Agreement  is  again  considered  good. 

Additional  calculations  were  performed  in  which  a  circular  bubble 
of  gas  was  collapsed  by  the  passage  of  a  shock  over  it.  (The  physical 
conditions  were,  of  course,  only  crudely  simulated,  since  two  materials 
were  not  possible  for  the  calculation  to  handle.)  The  results  showed 
no  particular  surprises,  and  since  no  quantitative  comparisons  could  be 
made,  the  results  are  not  shown  here. 
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V.  DIRECTION  FOR  FURTHER  DEVELOPMENT 


Besides  those  problems  which  have  already  been  brought  up  in  Parts 
I  to  IV,  there  are  several  others  which  come  to  mind. 

1.  The  method  should  be  adapted  to  cylindrical  coordinates.  Par¬ 
ticles  might  become  rings  around  the  central  axis  ( assuming  dependence  of 
the  flow  on  r  and  z  only),  and  the  appropriate  force  law  between  par¬ 
ticles  would  have  to  be  found.  If  a  ring  were  to  have  constant  mass, 
then  its  approach  to  the  axis  would  cause  difficulties;  thus  it  may  be 
necessary  for  ring  mass  to  be  time-dependent.  Boundary  conditions  will 
be  required. 

2.  Most  problems  of  interest  involve  several  materials  with  dif¬ 
ferent  state  equations.  A  method  must  be  evolved  for  treating  properly 
their  interaction. 

3.  The  effects  of  heat  conduction  and  true  viscosity  will  have  to 
be  considered;  it  is  not  expected  that  their  inclusion  will  be  difficult. 

Our  program  calls  for  an  exploration  of  these  and  other  matters  re¬ 
lative  to  eventual  development  of  a  useful  computing  method.  Comments 
and  suggestions  will  be  very  much  welcome  from  whomever  may  be  interested. 


APPENDIX  I.  CONFIGURATIONAL  INSTABILITY 


1  h 

Birkhoff  and  Lynch  have  remarked  that  a  rectangular  array  of  par¬ 
ticles  can  he  deformed  with  no  change  of  interparticle  distances,  suggest¬ 
ing,  therefore,  that  such  a  configuration  must  not  be  stable.  To  see  in 
more  detail  the  nature  of  configurational  instability,  we  consider  the 
simple  use  of  a  square  array  ...  .  . 


of  particles.  A  typical  par¬ 
ticle,  with  four  neighbors 
numbered  as  shown.,  will  have 
zero  net  force  on  it.  We  now 


suppose  that  all  five  particles  ...  .  .  . 

of  this  subset  be  displaced  by 

small  distances  without  change  of  internal  energy,  and  examine  the  re¬ 
sulting  force.  The  initial  coordinates  r*  and  (i  =  1 ,  2,  3,  L)  become 

J  1 

r*  +  I*  and  r\  +  e*.  and  the  resulting  force  is 
«3  d  i  i 


This  can  be  expanded  in  powers  of  €  and,  with  use  made  of  the  equilibrium 


conditions,  the  result  is 


— > 
F 


f  ^  \  f  \  f  j.  4  — > 

.  =  o  (4i>.  -  Z  t.  )  +  (  f  •  -  -2-)  4i*.  -  Z  (s. .) 

]  6xcA  J  i=i  ^  °  &v  l  J  i=i  ijo 


(S..)  •€, 
*• 


} 


where  prime  means  derivative  with  respect  to  argument;  6xq  is  the  unper¬ 
turbed  interparticle  distance;  and  fQ  =  f ( &Xq ) .  The  displacement  men¬ 
tioned  by  Birkhoff  and  Lynch  is  (to  lowest  order,  consistent  with  our  ex- 
pansions)  accomplished  by  taking  =  0  and  € ^  =  -€^  parallel 

to  a  line  from  #1  to  #3»  The  result  leaves  F.  =  0.  Various  other  dis- 

«J 

— ► 

placements  are  likewise  possible  with  F^  =  0. 

From  another  point  of  view,  we  may  examine  the  question  of  con¬ 
figurational  instability  as  follows.  As  a  result  of  the  displacements, 
if  the  forces  on  all  the  particles  are  directed  towards  returning  them 
to  their  initial  positions,  then  the  system  can  be  considered  stable.  We 
have  already  seen  that  such  a  tendency  to  return  cannot  in  general  result 
from  arbitrary  displacement.  There  is  one  special  case  of  interest,  how¬ 
ever,  that  is  easily  analyzed.  If  the  only  non-vanishing  displacement  is 
that  of  particle  then  the  resulting  force  on  is 

f  \ 


F.  =  2e.  (f  '  .+  jr2-  ) 
0  j  \  o  5x7 


and  the  condition  of  stability  is  thus 


f  '  +  <  0 

o  ox 

o 


This,  then,  is  a  necessary (but  not  sufficient)  condition  upon  the  force 
function  for  configurational  stability.  If,  for  example,  f(r)  is  pro- 
portional  to  r  b  then  the  stability  condition  is  £  >  1 .  For  a  polytropic 
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gas  in  adiabatic  motion,  this  will  be  satisfied,  while  in  isothermal 
motion  the  satisfaction  is  borderline. 


APPENDIX  II.  ANGULAR  MOMENTUM 


After  the  main  body  of  this  report  was  written,  it  was  discovered 

that  the  form  of  the  ef.  .  force  would  not  allow  angular  momentum  conser- 

ij 

vat ion,  and  a  review  of  the  situation  both  clarifies  the  angular  momen- 
turn  question  and  leads  to  an  alternative  expression  for  g  ..  No  quali- 

J- J 

tative  changes  would  result  in  the  discussions  of  this  report  from  the 
use  of  the  alternative  form  of  ~g,  . ,  since  the  two  forms  become  the  same 
in  one  dimension  and  differ  only  slightly  in  two. 

The  angular  momentum  of  particle  ifj  at  the  beginning  of  time  cycle 

#n  is 

r*  n  -»n  w  ->  n 
G  .  =  r .  X  m  .u . 

3  3  3  3 


so  that 


nn+1 

j 

3 

1  f-*n 

“2  (rJ 

->  n+1  s  v  /-+  n+1 

+  r .  )  X  mfy 

nx 

"ud  } 

/->  n+1 

+  (r3 

-  rj  >  X  2  "j(uj 

n+1 

+  u. 

J 

Now  if,  in  Eq.  (9)  of  the  text,  we  had  used  the  form 


(II-1) 


-»  n+1  -» n 

ll _ ~  J-  =  l({?n  +  £  n+1 )  (n-2) 

Bt  2V  j  3 
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for  advancing  the  particle  coordinates,  then  the  second  term  on  the  right 
of  (il-l)  would  vanish  and  angular  momentum  would  change  only  as  a  re¬ 
sult  of  the  moment  of  forces.  The  form  (II-2)  is  not  allowed,  however, 
because  of  its  undesirable  property  of  instability,  which  may  be  proved 
by  analysis  such  as  used  in  Part  II-B.  Thus  there  must  be  a  discrepancy 
in  angular-momentum  conservation,  with  cumulative  effect  proportional  to 
5t. 

We  may  also  examine  the  change  of  angular  momentum  of  any  subset  of 
particles,  and  find  for  conservation  (to  lowest  order  in  5t)  that  the 
interparticle  forces  must  lie  along  the  lines  of  centers.  This  leads  to 
the  alternative  form  for  the  dissipative  force 
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